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Abstract 

o 

■ The most important goal of studying an extensive air shower is to find the en- 

ergy, mass and arrival direction of its primary cosmic ray. In order to find these 
parameters, the core position and arrival direction of the shower should be deter- 
mined. In this paper, a new method for finding core location has been introduced 
that utilizes trigger time information of particle detectors. We have also developed 
a simple technique to reconstruct the arrival direction. Our method is not based 
upon density-sensitive detectors which are sensitive to the number of crossing par- 
ticles and is also independent of lateral distribution models. This model has been 
developed and examined by the analysis of simulated shower events generated by 
the CORSIKA package. 
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1 1 Introduction 



2 An Extensive Air Shower (EAS) is a large number of secondary particles orig- 

3 inating from a high energy primary cosmic particle. A lot of efforts have been 

4 made to understand the structure and development of EASs. For an accu- 

5 rate investigation of the EAS events, we first have to know their direction 

6 and core position. The better the measurement of these parameters, the more 

7 precise the exploration of extensive air shower structures. Plane Wave Front 

8 Approximation (PWFA) is the simplest way to find the direction of EAS, and 
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1 spherical front approximation [1] is a complementary approach to obtain more 

2 accurate results. The common method for finding core locations of extensive 

3 air showers is to fit a Lateral Density Function (LDF) to the density of sec- 

4 ondary particles of the shower. Greisen function [2] and its modified form (e.g. 

5 [3]) are usually used as LDF. In LDF methods, one uses the particle density 

6 information without considering their arrival time information. 

7 Some attempts have been made to use mean arrival time and disk thickness 

8 as a measure of the shower core distance [4]. Measurements of the EAS disk 

9 structure were tried by Bassi, Clark, and Rossi in 1953 [5] and continued later 

10 [6]. However, it has recently been shown that both mean arrival time and EAS 

11 front thickness in individual showers fluctuate strongly and cannot be a good 

12 measure of the distance from the EAS axis [7]. 

13 Though, due to the strong fluctuations, arrival time cannot be used as a mea- 

14 sure of the distance from the core, it can be used to provide a statistical 

15 analysis for some of the EAS features. In this paper, we propose a model inde- 

16 pendent approach for finding core location and arrival direction of EASs which 

17 uses arrival time information, avoids using LDF and can be used for arrays 

18 lacking Density Sensitive Detectors (DSDs), which are sensitive to density or 

19 to the number of crossing particles. In order to examine the capability of the 

20 method, we used the CORSIKA [8] package. 



21 2 Physical Principles 

22 An important feature of EASs is their spherical front, which is approximately 

23 a spherical cap. Thus the particles in the core region of a Vertical EAS (VEAS) 

24 reach to the ground level sooner than in other regions. This feature can be 

25 demonstrated by considering the fact that if we randomly select any two sec- 

26 ondary particles of a VEAS, the first particle reaching ground level on average 

27 is closer to the core location. Fig. 1 shows that by increasing the distance be- 

28 tween two particles, the average distance of the first arriving particle from the 

29 core increases slowly, while the average distance of the second particle from 

30 the core increases rather rapidly. 

31 Another feature of EASs is that the particle density in the core neighborhood 

32 is larger than in other regions. The smaller the distance between two particles, 

33 the smaller their distance to the core on average, as also follows from Fig. 1. 

34 Error bands also have been depicted in order to show the fluctuations. They 

35 have been derived by calculating the RMS for positive and negative errors 

36 separately. As can be seen the error bands are asymmetric. Fig. 2 shows the 

37 reason of this asymmetry. 

38 The peak structures in Fig. 2 show that the first arriving particle to the ground 

39 level is more likely to be closer to the core location. For the small separations 

40 of particles, both probability density functions have a unique pronounced peak 
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1 close to the core region. So, when the separation of two particles is small, it is 

2 quite probable that both particles to be in approximately equal distance from 

3 the core in opposite sides. For larger separations, both probability densities 

4 have two peaks, one of which is higher than the other. For the first arriving 

5 particle, the higher maximum is near the core region, but the other one has 

6 the same distance from the core as the distance between two particles. For the 

7 second particle the situation is vice versa. So, when their distances increase, 

8 although there is a small probability that the first arriving particle will be 

9 found farther from the core than the second one, it is more probable that the 

10 first arriving particle to be much closer to the core region than the second 

11 arriving particle. 

12 The density of secondary particles quickly decreases by increasing the distance 

13 from the core in contrast to the slow increase of the average arrival time of 

14 particles shown in Fig. 3. A rough investigation shows that random fluctuation 

15 of the density of the particles is much less than their arrival time fluctuations, 

16 especially when they are far from the core. Fig. 4 shows that the fluctuation 

17 of lateral density of the particles decreases with the core distance rapidly in 

18 contrast to the arrival time fluctuation which increases with the core distance 

19 (Fig. 3). So if we want to choose between particle density information and 

20 arrival time information as a measure of the core location, the particle density 

21 information is a more decisive factor. 



22 3 A Method for Finding Core Location 

23 A simple approximate method to find the core position of an EAS is to cal- 

24 culate the Center of Gravity (CG) of the Triggered Detectors (TDs) which is 

25 a relatively good approximation for those EASs whose cores are very close to 

26 the center of the array. We have introduced a procedure which is capable to 

27 increase the precision of CG for finding the approximate location of the core 

28 even when the core position is close to the border of the array. 

29 At first, we assume that every detector can just detect the first crossing parti- 

30 cle and is unable to detect any other particle during an event. For simplicity, 

31 we investigated VEAS events at first and then generalized the results to the 

32 inclined EASs. 

33 Assume that we have positions and trigger times of all TDs of an array during 

34 a VEAS event. TDs are indexed based on their trigger times, starting with 

35 % — 1 for the first TD. Fig. 3 suggests that if % < j < k then (r^) < (r,) < (r k ), 

36 where {r\) is the average distance of the ith TD to the shower core. Therefore, 

37 by taking into account Fig. 1, this result can be achieved: (dij) < {dik) (where 

38 d^ = \ fi — fj\). If we want to find the nearest detector to the core, selecting 

39 the ith detector seems to be logical, but because of timing fluctuations, a bet- 

40 ter choice will be obtained by the following procedure: At first, we find the 
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1 minimum value of dij, do- and djk- If djk is the smallest value then the jth 

2 detector is most likely nearest to the core, since j < k. 

3 3.1 SIMEFIC Algorithm 

4 Based on the principles explained above, we developed the SIMEFIC (Sieving 

5 MEthod for Finding Core) algorithm for eliminating detectors far from the 

6 core. Let us assume that there are N TDs in an array during a VEAS. We 

7 form a matrix D^xn, whose elements are d^s. In view of the fact that D is a 

8 symmetric matrix, we just consider the upper triangle of the matrix without 

9 principal diagonal elements, which are zero. 
10 
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11 Then we can find the smallest element (d<) under the following two conditions: 

12 If d^ = dki and % < k then d< = d^ 

13 If d^ = dik and j < k then d< = d^ 

14 Now we select the ith TD as the first detector of our near-core list (because of 

15 its trigger time). Next, in the ith row, we find the biggest element (d>) under 

16 the following condition: 

17 If all dijS (j — % + 1, N) are different, the biggest element, e.g. dik, is d>, 

18 and if d^ and d^ are both the biggest elements (d^ = dik), and j < k then 

19 d > = d ik . 

20 We now eliminate the fcth TD as an off-core detector and remove the ith and 

21 kth. rows and columns of the matrix D . By repeating this procedure for the 

22 reduced matrix (Djv-2xiv-2), we will reach a position in which half of the TDs 

23 are retained and half of them are eliminated. Now, it is expected that the CG 

24 of the retained TDs is a good measure of the core location. 

25 3.1.1 Inclined showers 

26 Up to now, our discussion was limited to the VEASs. To generalize this method 

27 to inclined EASs, the detector coordinates, and also the trigger times need to 
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1 be transformed into the coordinate system of the shower. The coordinates of 

2 the detectors on the ground level are Xi, yi, and Zi = 0. Now, these coordinates 

3 are transformed to a new coordinate system whose x and y axes are perpen- 

4 dicular to the arrival direction. To perform this, we select an arbitrary point 

5 (e.g. the CG of all detectors) as the origin of coordinate system. Then the 

6 coordinates of the detectors are first rotated counterclockwise around the z 

7 axis by the angle and next the new coordinates are rotated clockwise around 

8 the y by the angle 9. The z components in the final coordinate system are 

9 z i = —Xi sin 6 cos <t> — Vi sin 6 sin 0. Accordingly the trigger times of the detec- 

10 tors are transformed to t'- = U + z- /c where c is the speed of light. Now the 

11 inclined showers are treated as the vertical ones. 



12 3.1.2 Arrays of DSDs 

13 For arrays with DSDs the algorithm is extended as follows. Another matrix, 

14 N, is formed, the elements of which are riij = rii + rij (rii is the number of the 

15 detected particles by ith detector) in correspondence with dij elements of the 

16 matrix D. Now, we first select a pair of detectors with maximum (n max ) 

17 under the following two conditions: 

18 If riij = riki with i < k, and if > dki, then n max = riki and if d^ < dki, 

19 then n max = n^. 

20 If riij = riik with j < k, and if d^ > da-, then n max = and if d^ < dik, 

21 then n max = n^. 

22 Assume that riy has been selected as n max and rii < n j, then j is the near 

23 core detector (if rii > n-j, «th detector will be chosen). If rii = n-j, the trigger 

24 times of the two detectors % and j will be the determining factor for choosing 

25 between % and j as the near core detector. Other stages of the algorithm are 

26 the same as before. 



27 3.2 Reconstruction of a Shower Geometry 



28 As the precision of final results of SIMEFIC algorithm is tied to the precision 

29 of arrival direction measurement, a method has been introduced here that can 

30 be used for more precise arrival direction reconstruction. 

31 Assume that we have {Pj} set of the locations and trigger times of all TDs at 

32 the beginning of the calculation. By using PWFA and SIMEFIC algorithm on 

33 the set {Pi}, a subset of it, {P-}, will be found whose members are half of the 

34 first set. Now, PWFA is used in shower direction reconstruction of the {Pi} 

35 set from which the arrival direction is calculated with more precision. The 

36 reason for higher precision, of this technique in comparison with the PWFA 

37 used for the set {Pi}, is that arrival times have less fluctuation in near core 
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1 region than those in far regions from the core. Furthermore, the front of an 

2 EAS is smoother in near core region than that in far regions. In view of these 

3 facts, SIMEFIC method presents a better approximation for arrival direction 

4 in comparison with the PWFA, which uses data of all TDs. Now, we propose 

5 a refinement to the SIMEFIC method by repetition of the above procedure: 

6 (1) Reconstruct arrival direction by using {PJ set with PWFA, (0,0). 

7 (2) Use SIMEFIC algorithm for the original set, {p} (not on set {Pj } which 

8 has N/2 members) and arrival direction (#',0'), to obtain another set 

9 W>- 

10 (3) Replace {P t } by {P i } and repeat the procedure. 

11 This procedure is repeated several times using the original set to refine the 

12 arrival direction. If two successive runs yield approximately the same results, 

13 the process will be terminated. In our simulation, repeating for 3 times was 

14 enough. 



15 4 Simulations with CORSIKA 



16 The geographical coordinates of the prototype of Alborz observatory loca- 

17 tion (in Tehran, 35°N, 51°E and 1200 m above sea level) have been imposed 

18 in CORSIKA EAS simulations. Geomagnetic field components of Tehran are 

19 B x = 28.1/iT and B z = 38.4yuT. For the simulation of low energy hadronic 

20 interactions, GHEISHA [13] package, and for high energy cases QGSJET01 

21 package have been used. Zenith angles of primary particles were chosen be- 

22 tween 0° and 60°. The compositions of primary particles have been 90% pro- 

23 tons and 10% helium nuclei. This ratio was assumed to be constant over entire 

24 energy range. The energy of the primary particles ranges from 100 TeV to 5 

25 PeV. Other parameters are CORSIKA default values (e.g. default spectrum 

26 index of -2.7). 

27 The assumed array is a square with an area of 200 x 200 m 2 , composed of 

28 1 x 1 m 2 detectors on a square lattice with a 5 m lattice constant (total of 

29 41 x 41 detectors). Ground arrays are commonly made up of scintillation de- 

30 tectors which are not sensitive to the position of the crossing particles and are 

31 not able to exactly measure how many particles have passed through them. 

32 Therefore, the following assumptions have been applied for the simulation. 

33 When the first arriving particle passes through a detector, its arrival time is 

34 considered as the trigger time of that detector. Obviously, the coordinates of 

35 this first crossing particle must be within ±0.5 m of the center of the detector 

36 in order to consider the detector as a TD. The next arriving particles which 

37 pass through this hypothetical detector are not taken into account. 

38 The CG of TDs for the array whose detectors are DSDs will be considered as 
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1 follows: 

x CG = l/N^mxi, 

i 

yea = I/N^mvi 

i 

2 where rij is the pulse height (the number of crossing particles) of each detector 

3 and iV = ^ n^. For non DSDs, we set n« = 1 for all TDs. 

i 

4 If the coordinates of the EAS true core are denoted by (x tc ,ytc), the distance 

5 between the true core and the CG will be: 

r = \J {x' CG - x't c ) 2 + {y' CG - y'tc) 2 

6 where x",y" are coordinates in the rotated coordinate system introduced in 

7 sec. 3.1.1. 
8 



9 4-1 Results 

10 Fig. 5 shows an EAS which its true core is on point (—80,-80). Detectors 

11 accepted by SIMEFIC method are shown by the bold circles and the omitted 

12 ones by the empty circles. It is clear that even when the true core is near the 

13 edge of the array, this algorithm has approximately chosen proper detectors, 

14 while, some of the detectors that are far from the core but near to each other 

15 have not been chosen. This is an important effect of using the timing infor- 

16 mation. This method is precise up to the point that almost half of the TDs 

17 are symmetrically spread around the core. 

18 In Fig. 6 and Fig. 7, we consider the true core on the diagonal line on the 

19 points (10i, 10i) (i = 1,...,10), with the center of the array at point (0,0). 

20 Trigger condition for the EASs, whose specifications have been introduced at 

21 the first part of the current section, was triggering at least 68 detectors out of 

22 1681 detectors by secondary particles. 3000 of the accepted EASs have been 

23 averaged for each data point shown in Fig. 6 and Fig. 7 and since the size 

24 of error bars for each data point is less than the size of the symbols used for 

25 them, the error bars are omitted in these figures. 

26 Fig. 6 shows the results of using SIMEFIC method and the CG of all TDs for 

27 finding core locations. It is clear that the CG of all TDs is only precise for 

28 those EASs whose core are near the central part of the array and, by increas- 

29 ing the distance of EAS core from the center of the array, the error increases 

30 gradually. But, in SIMEFIC method, the precision, even up to 50% of the 

31 length of the array, is approximately the same. As we expect, the results for 

32 an array with DSDs are better than those for the array without DSDs. 
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1 Fig. 7 shows the angle between the EAS primary direction provided by COR- 

2 SIKA and the direction which has been found by the SIMEFIC method and 

3 also by PWFA (applied on all TDs) versus the distance of true core from 

4 array center. It is clear that SIMEFIC algorithm has significantly improved 

5 PWFA. Furthermore, Fig. 7 shows that the results obtained with non DSDs 

6 are slightly more precise than those of DSDs. We think that the reason behind 

7 this unexpected result is that for finding direction of an EAS, it is better to 

8 give all of the detectors the same weight and choose them with the same prob- 

9 ability around the true core location. When we discard density information, 

10 all of the TDs have the same probability to being chosen, so the algorithm 

11 will select TDs symmetrically around the core. Due to the stochastic nature 

12 of secondary particles of the shower, the symmetric selection will be rather 

13 destroyed with DSDs and those TDs which have more detected particles have 

14 a more chance for selection. 



15 5 Conclusion 



16 In this investigation, we have developed a simple method for finding the posi- 

17 tion of an EAS core, and applied PWFA to an optimized data set for finding 

18 the EAS arrival direction. In this method, we have used the information of 

19 positions of the TDs on the ground level as well as the trigger times. In this 

20 method, distances between all pairs of TDs are measured and then TDs with 

21 minimum separation and smallest trigger time are chosen, while TDs with the 

22 maximum separation and largest trigger times are omitted. Finally, the CG 

23 of the chosen TDs is used for finding core location of EAS. 

24 The essence and scheme of the method envisioned by examining vertical show- 

25 ers have been also formulated for inclined showers. An operational algorithm 

26 were developed and tested over 3 x 10 5 simulated EAS events generated by 

27 CORSIKA package. The proposed analysis technique is adequate for simple 

28 EAS arrays without DSDs, though it has been generalized for arrays with 

29 DSDs, and its results for finding core location are more accurate with DSDs. 

30 The precision of finding EAS core location is improved to about the array's 

31 lattice constant for an EAS whose core falls within the central region of the 

32 array and to about 4 times lattice constant for those falling close to the edge 

33 of the array. An improvement of about 58% has been reached in comparison 

34 with the CG of all TDs for the above mentioned geometry and configuration 

35 of the array. Furthermore, by using the PWFA for the TDs selected by this 

36 method, the angular resolution of the primary arrival direction is improved 

37 significantly compared to the case of using PWFA on all TDs of the array 

38 during an event. The angular error ranging from 1.2° to 4.2° for the case of 

39 using all TDs is reduced to values ranging from 0.4° to 1.5° corresponding to 

40 showers falling near the center or on the edge of the array. Again, this is on 
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1 average an improvement of about 89% in comparison with the simple PWFA. 

2 It should be noted that these improvements belong to the results of our sim- 

3 ulations for the assumed array of 1681 normal detectors. Obviously, the im- 

4 provement depends on the size of EAS array, its spacing, and its detectors 

5 type. Further improvements are expected by optimally combining the infor- 

6 mation contained in the two matrices, D and N mentioned in sec. 3.1. In our 

7 future attempts we will investigate simulated data to find the optimal method 

8 of combining the information of these two matrices (D and N). The method 

9 for finding core location could be used as an improved first guess in order to 
10 seed the common fitting methods. 

11 
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Fig. 1. Average distance of the first arriving particle (solid curve) and the second 
arriving particle (dash-dot curve) from the core versus their distance from each 
other for all secondary charged particles of the showers. The related error bands 
are also depicted. The results have been reached by additional 10,000 vertical single 
energy (100 TeV) showers generated by CORSIKA. The low energy hadronic model 
is FLUKA [9,10] and the high energy hadronic model is QGSJETII [11,12]. 
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Fig. 2. Probability density of distances of two typical particles of those EASs which 
were used in Fig. 1. Three groups of particle couples were selected with separations 
of 20 m ± 0.5 m, 50 m ± 0.5 m and 100 m ± 0.5 m. Solid lines belong to the first 
arriving particles and dashed lines belong to the second arriving particles. In all 
three cases the first arriving particle is more probable to be closer to the core. 
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Fig. 3. Mean arrival times of particles versus the distance from the core for the same 
showers which were used in Fig. 1. Dashed lines are error bands. 
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Fig. 5. The TDs which are accepted by SIMEFIC method are shown by the bold 
circles and the omitted ones by the empty circles. The position of the true core and 
the center of the array have been shown. The zenith and azimuth angles of this 
shower are 29.5° and 329.5° respectively, and its energy is 157.5 TeV. 
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Fig. 6. The comparison of the CG of all TDs and those which have been selected 
by SIMEFIC. The size of error bars for each data point is less than the size of the 
symbols used for them. 
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Fig. 7. Mean value of angle between the EAS primary direction and the direction 
which has been found by SIMEFIC algorithm assuming the array with non DSDs 
(solid curve), SIMEFIC algorithm assuming the array with DSDs (dashed curve) 
and by PWFA applied on all TDs (dotted curve) versus the distance of the true 
core from array center. 
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